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METHOD THEREFOR 



Field of the Invention 

The present invention relates to color imaging. 
More particularly, the present invention relates to mapping 
colors between color imaging systems . 

Background of the In vention 

Color reproduction processes typically involve 
using color imaging systems to produce colors on various 
media. These color imaging systems may be used to duplicate 
a color image from one medium to another medium, e.g., from 
one printed copy to another or from a display screen to a 
printed copy. Color reproduction processes are used in 
various 'application environments, for example, color 
proofing applications. 

Some color reproduction processes use approaches 
known as color management systems (CMSs) to characterize 
various color imaging systems and to transform color data 
between the color imaging systems. Characterizing color 



imaging systems typically involves calculating color 
response functions using color coordinate systems known as 
color spaces. One commonly-used color space is Commission 
Internationale de l'Eclairage L*a*b* (CIELAB) space. CMSs 
5 attempt to reproduce an original color image on a color 
imaging system so as to preserve the appearance of colors 
between the original and the reproduction within the 
limitations of the color imaging system of the reproduction 
process . 

10 Various CMS approaches have been proposed to 

achieve accurate color reproduction. Many of these 
approaches involve producing color samples using an output 
or display device and measuring the color values of the 
samples using an input device. Such approaches correlate 
15 the output colors with the measured color values. This 
correlation is performed using, for example, forward and 
reverse transforms between device-specific color spaces and 
a device -independent color space. These transformation 
techniques are often supplemented by interpolation between 
2 0 entries in a multidimensional lookup table. These 

techniques exhibit inaccurate color conversion between 
similar devices, potentially resulting in undesirable 

2 



contamination of colors. Furthermore, accurate color 
conversion of dark colors has often been particularly 
difficult because of inadequate processing of black channel 
data in many applications. 
5 CMSs often perform gamut mapping to correlate the 

range or gamut of colors that can be realized by a device 
with regions of a color space. Because many devices are 
incapable of realizing the complete range of colors in a 
color space, gamut mapping typically involves compressing or 
10 scaling regions of the color space. The device can then 
approximate colors outside its gamut using the compressed 
regions of the color space. For many CMSs, gamut mapping is 
potentially inconsistent under certain circumstances, such 
as when using profiles generated by software from different 
4 15 vendors. In addition, many CMSs exhibit inconsistencies 

when performing forward and reverse transformations between 
imaging systems. For example, color shifting often occurs 
with repeated forward and reverse transformations. 

Many CMS techniques exhibit other limitations in 
20 addition to the lack of accuracy in converting colors. For 
example, many CMS techniques are relatively inflexible with 
respect to changes in illumination and observer conditions, 
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gamut mapping, and choice of color space. Certain 
techniques lack forward compatibility with future color 
standards . 

5 Summary of the Invention 

According to one embodiment, the present invention 
is directed to a color mapping method for use in 
transforming colors between color imaging systems. The 
color mapping method includes using forward transformation 
10 profiles characterizing the color imaging systems to 

generate respective sets of device -independent color values 
for the color imaging systems. Color conversions are 
calculated by recursively reducing differences between the 
sets of device- independent color values. This difference 
• 15 reduction is also optionally performed on black channel 

information to obtain a mapping of black channels between 
the color imaging systems . A color map describing a 
relationship between the color imaging systems is 
constructed as using the predicted color conversions. This 
20 method may be performed by a color mapping arrangement or a 
computer-executable program stored on a data storage medium. 



According to another embodiment of the present 
invention, color mapping between source and destination 
color imaging systems is accomplished by using profiles that 
characterize the color imaging systems to generate device- 
5 independent color values for the source color imaging system 
and to convert to device -dependent values of the destination 
color imaging system by performing a color conversion using 
the profiles. The device -independent color values have a 
same dimensionality as the corresponding color imaging 

10 systems. The color conversion can be used to improve its 
own accuracy relative to a quality threshold. The color 
conversion is used to define a color map for transforming 
colors between the color imaging systems . 

Another embodiment of the present invention is 

15 directed to a color mapping arrangement for use in 

transforming colors between imaging systems. A computer 
arrangement uses forward transformation profiles that 
characterize the color imaging systems to generate 
respective sets of device- independent color values for the 

2 0 color imaging systems. The computer arrangement also 
calculates color conversions by recursively reducing 
differences between the sets of device -independent color 
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values. The computer arrangement uses the color conversions 
to construct a color map describing a relationship between 
the color imaging systems using the color conversions. A 
memory stores the color map. 
5 The above summary of the invention is not intended 

to describe each disclosed embodiment of the present 
invention. This is the purpose of the figures and of the 
detailed description that follows. 



10 Brief Description of the Drawings 

Other aspects and advantages of the present 
invention will become apparent upon reading the following 
detailed description and upon reference to the drawings in 
which: 

* 15 FIG. 1 is a block diagram illustrating an example 

color mapping system, according to an embodiment of the 
present invention; 

FIG. 2 is a block diagram illustrating an example 
arrangement implementing part of the color mapping system of 
20 FIG. 1, according to an embodiment of the present invention; 

FIG. 3 is a block diagram illustrating another 
example arrangement implementing part of the color mapping 
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system of FIG. 1, according to an embodiment of the present 
invention; 

FIG. 4 is a block diagram illustrating yet another 
example arrangement implementing part of the color mapping 
5 system of FIG. 1, according to an embodiment of the present 
invention; 

FIG. 5 is a block diagram illustrating still 
another example arrangement implementing part of the color 
mapping system of FIG. 1, according to an embodiment of the 
10 present invention; and 

FIG. 6 is a flow chart illustrating an example 
color mapping method, according to another embodiment of the 
present invention . 

While the invention is amenable to various 
* 15 modifications and alternative forms, specifics thereof have 
been shown by way of example in the drawings and will be 
described in detail. It should be understood, however, that 
the intention is not to limit the invention to the 
particular embodiments described. On the contrary, the 
20 intention is to cover all modifications, equivalents, and 
alternatives falling within the spirit and scope of the 
invention as defined by the appended claims. 
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Detailed Descrip tion of the Various Embodiments 

The present invention is believed to be applicable 
to a variety of systems and arrangements that characterize 
5 color imaging systems. The invention has been found to be 
particularly advantageous for transforming colors between 
different color imaging systems. An appreciation of various 
aspects of the invention is best gained through a discussion 
of these particular application examples. 

10 According to one aspect of the present invention, 

a color mapping technique may be applied to a variety of 
color imaging systems to generate a color map that can be 
used to transform the color response of one color imaging 
system, referred to as a source color imaging system, to 
* 15 match the color response of another color imaging system, 
referred to as a destination color imaging system. The 
color mapping technique projects color coordinates in the 
color space used by the source color imaging system into, 
for example, a device -independent color space. Optimal 

20 color coordinates in the color space used by the destination 
color imaging system are determined that realize a 
relatively close match between the projections into the 
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device- independent color space of the color coordinates in 
the color space used by the source color imaging system and 
the optimal color coordinates. The color mapping technique 
then generates a color map based on the optimal color 
5 coordinates for a number of color coordinates in the color 
space used by the source color imaging system. 

FIG. 1 illustrates an example system 100 according 
to the present invention configured to transform colors 
between imaging systems. The system 100 includes an 

10 appropriately-programmed computer arrangement 102. The 

computer arrangement 102 may be implemented using any of a 
variety of conventional resources, for example, a personal 
computer and CD-ROM based software. Other computer-based 
designs may be used as well. For example, the computer 

15 arrangement 102 may be implemented using a microprocessor 

that acts as a read-only memory (ROM) into which a software 
application program is loaded. The software application 
program may be incorporated, for example, in a color- 
management software package . 

20 A color mapping system 104 includes a color 

management system 106. The color management system 106 
receives a source device profile 108 and a destination 



device profile 110. These device profiles describe mappings 
from device-dependent color coordinate systems used by 
respective color imaging systems to device- independent color 
coordinate systems . 
5 The color management system 106 processes the 

source device profile 108 and the destination device profile 
110 to generate a color map 114. The color map 114 
describes a relationship between the color imaging systems 
used by the source and destination devices. A memory 116 

10 stores the color map 114. Subsequently, the color 

management system 106 uses the color map 114 to transform a 
set of source coordinates 118 in a device -dependent source 
device color space into a set of destination coordinates 12 0 
in a device -dependent destination device color space. 

15 FIG. 2 illustrates an example color management 

system 20 0 for transforming colors between imaging systems 
according to the present invention. A source device profile 
interpreter 202 receives a source device profile 206. The 
source device profile 206 is used to map coordinates in the 

2 0 source device color space to a some form of color data, such 
as spectral or XYZ tristimulus values. For example, if the 
source device is a halftone color printer, the source device 



profile 2 06 may map CMYK color values to a XYZ color space. 
The source device profile interpreter 2 02 interprets the 
source device profile 206 and converts coordinates in the 
source device color space to a device -independent color 
5 space known as a profile connecting space (PCS) . The PCS is 
used for converting the coordinates in the source device 
color space to the destination device color space. The PCS 
may be, for example, the CIELAB color space. Another 
example PCS is described in copending U.S. Patent 
10 Application, entitled ''Characterization of Color Imaging 

Systems" (Christopher Edge et al.), assigned to the instant 
assignee, filed on June 27, 1997, and incorporated herein by 
reference . 

A destination device profile interpreter 208 
15 receives a destination device profile 210. The destination 
device profile 210 is used to map color coordinates in a 
destination device color space used by a destination device 
212 to some form of color data, such as spectral or XYZ 
tristimulus values. For example, if the destination device 
20 212 is a cathode ray tube (CRT) monitor, the destination 

device profile 210 may map color coordinates in a red-green- 
blue (RGB) color space to XYZ tristimulus values. The 
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destination device profile interpreter 208 interprets the 
destination device profile 210 and converts color 
coordinates in the destination device color space to the 
PCS. 

5 The source and destination device profile 

interpreters 202 and 208 may be implemented using any of a 
variety of hardware and software arrangements and are 
configurable for a variety of application environments. For 
example, if the source and destination device profiles 2 06 
10 and 210 are International Color Consortium (ICC) device 
profiles, the source and destination device profile 
interpreters 202 and 208 are optionally configured to 
include white- and black-point parameters to account for 
color variations between media and colorants used by 
-15 different color display devices. The source and destination 
device profile interpreters 202 and 208 can also be 
configured to include pleasing color corrections, such as L* 
rescaling and a*b* hue adjustments. Alternatively, the 
pleasing color corrections can be incorporated into the 
20 color transformer 214. In certain other application 

environments, the source and destination device profile 
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interpreters 202 and 208 are further configurable to 
include, for example, illuminant and observer functions. 

The device profile interpreters 202 and 2 08 can be 
configured using any of a variety of approaches. For 
5 example, plug- in software modules can be used to configure 
the device profile interpreters 202 and 208. Using plug-in 
software modules obviates the need to use new versions of 
the color management system 200 or of the device profiles 
2 06 and 210 when adding, for example, a newly defined color 

10 space, a custom illuminant, such as fluorescent light, or a 
new gamut mapping technique. These options can be selected, 
e.g., using a setup window at the operating system level. 
For example, if the operating system is Apple OS version 
7.5, these options can be selected using a control panel 

15 interface . 

If the device characterization is non-spectral , 
the color management system 200 can use the original 
spectral data that is saved with the profile to reconstruct 
the device profiles according to various conditions, such as 
20 illuminant functions and color space choices. For example, 
if one uses an RGB regression to convert scanner RGB values 
into color space values for a particular combination of 
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color space and illuminant and observer conditions based on 
a set of spectral data, the regression for a new set of 
conditions can be generated based on the same spectral data. 
Accordingly, the device profiles 206 and 210 can be used to 
5 calculate color values for a variety of conditions and color 
appearance models * 

A color transformer 214 obtains PCS color 
coordinates from the source and destination device profile 
interpreters 202 and 208. The color transformer 214 uses 

10 these color coordinates to develop a color map 216 that 

expresses a relationship between the color spaces used by 
the source and destination devices 204 and 212. To generate 
the color map 216, the color transformer 214 may use any of 
a variety of gamut mapping techniques. One such technique 

15 that has been found to yield particularly accurate results 
involves reducing the color error between the source and 
destination devices. The color error is defined, for 
example, by Euclidean distances in the PCS or by weighted 
sum square errors in a color space that is polar in the 

2 0 chromatic dimensions of the PCS. Defining the color error 
using weighted sum square errors results in a mapping 
between color imaging systems that accurately maintains 



colors in reproduced images. By using error reduction 
techniques, the color transformer 214 avoids generating 
significant cumulative error in performing multiple forward 
and reverse transformations between color spaces. 
5 The color transformer 214 is implemented using, 

for example, a software program, and can be configured for a 
variety of applications. For example, the color transformer 
214 can be configured to perform a 100% black point scaling 
for mapping a printed color image to a monitor display of 

10 the image. On the other hand, because newsprint has a 

relatively weak black point attributable to its ink density 
and light- transmitting properties, the color transformer 214 
can be configured to perform, for example, a 50% black point 
scaling when mapping a color image printed on newsprint to 

15 the Matchprint™ color imaging system. The color 

transformer 214 is also configurable to use, for example, 
illuminant and observer functions, which the color 
transformer 214 provides to the source and destination 
device profile interpreters 202 and 208. The color 

2 0 management system 2 00 receives user preferences from an 
input 218 to determine how to configure the color 
transformer 214. 



After developing the color map 216, the color 
transformer 214 can be used to transform colors between the 
source and destination devices 204 and 212. The color 
transformer 214 receives color coordinates from the source 
5 device 204 and transforms them using the color map 216. 

This transformation produces a set of color coordinates in 
the destination device color space. The destination color 
imaging system then reproduces the color on the destination 
device 212 using these color coordinates. 

10 FIG. 3 illustrates an example device profile 

interpreter 3 00 implementing part of the color management 
system 200 of FIG. 2. The device profile interpreter 300 
uses a device profile 3 02 to convert device coordinates 
received at an input 3 04 to PCS color coordinates, which the 
* 15 device profile interpreter 300 provides at an output 306. 
The device profile 3 02 describes the relationship between 
the device coordinates and some form of color data. 
Additionally, the device profile 302 optionally stores the 
raw spectral data used to construct the device profile 302. 

20 The raw spectral data allows subsequent construction of more 
accurate device profiles 302, e.g., if ICC specifications 
change. This updating can be performed automatically, for 
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example, upon detecting that some component of the device 
profile 3 02 is out of date. Updates can also be performed 
periodically based on a schedule. To update the device 
profile 302, a new profile can be generated using the 
5 spectral data. Alternatively, error reduction, such as a 

one -dimensional correction, can be performed on each channel 
in the original look-up table for constructing the new 
profile. This correction can be applied as a separate set 
of one-dimensional tables or applied directly to the 

10 analytical model or multidimensional look-up table. For 

additional information concerning an example error reduction 
procedure that can be used in constructing a new profile, 
reference can be made to U.S. Patent Application Ser. No. 
08/431,614, entitled "Apparatus and Method for Recalibrating 
• 15 a Multi-Color Imaging System," assigned to the instant 
assignee and incorporated herein by reference. 

A device profile processor 308 receives the device 
coordinates from the input 304 and the device profile 302. 
The device profile 3 02 may be, for example, an ICC profile. 

20 If the device profile 302 exists in this format, the device 
profile processor receives the forward portion of the 
profile, i.e., the portion used for converting device 
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coordinates to PCS color coordinates. Alternatively, the 
device profile 302 can be stored in another format. The 
device profile processor 3 08 processes the device 
coordinates using the device profile 3 02 and outputs certain 
5 data based on the device profile 302. For example, if the 
device profile 302 is an ICC profile, the device profile 
processor 3 08 outputs XYZ tristimulus values for a 
particular set of observer conditions {e.g., illuminant and 
observer functions) . If the device profile 302 is based on 

10 spectral data, the device profile processor 3 08 outputs 
spectral data. The device profile processor 3 08 can be 
configured for a variety of applications. For example, a 
user can select between absolute and relative colorimetrics 
and can configure observer, e.g., illuminant, conditions. 

15 A PCS processor 310 receives the data output from 

the device profile processor 3 08 and a set of PCS parameters 
from an input 312. The PCS parameters may include, for 
example, XYZ tristimulus values for the media white, the 
illuminant white, and the black point, as well as black- 

2 0 point scaling from a perfect black to the media black. The 
PCS processor 310 generates the PCS values as a function of 
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the data received from the device profile processor 3 08 and 
the PCS parameters. 

FIG. 4 illustrates an example color transformer 
400 implementing part of the color management system 200 of 
FIG. 2. A device link generator 4 02 receives as input at 
least one source profile and one destination profile. While 
FIG. 4 illustrates the device link generator 4 02 receiving a 
source device profile from a source profile interpreter 404 
and a destination device profile from a destination profile 
interpreter 4 06, it should be understood that the device 
link generator 4 02 may also receive one or more device 
profiles that are intermediate between the source and 
destination device profiles. For example, a device profile 
characterizing an RGB monitor can be intermediate between a 
source device profile characterizing an RGB scanner and a 
destination device profile characterizing a CMYK printer. 
The source and destination device profiles are forward 
transforms and optionally include configurable observer 
conditions and PCS parameters. The device link generator 
402 also receives a series of PCS parameters 408 to improve 
linking of different device types (e.g., CRT monitors and 



printers) . The gamut mapping parameters 410 improve mapping 
of out of gamut colors between device types . 

The device link generator 4 02 generates a color 
map or device profile link 412 that maps colors between two 
5 devices, e.g., from an RGB device to a CMYK device or 

between two CMYK devices. The device profile link 412 is, 
for example, a mathematical expression or a look-up table. 
The color transformer 400 optionally stores the device 
profile link 412 in a memory, such as a random access memory 
10 (RAM) , or saves it as a file for multiple transformations 
between the source and destination device color spaces. 

A device link calculator 414 receives source 
device coordinates from an input 416 and processes them 
using the device profile link 412. The device link 
• 15 calculator 414 uses a single forward calculation to 
transform the source device coordinates to a set of 
destination device coordinates for presentation at an output 
418. Because the device link calculator 414 uses a single 
forward calculation, interpolation is relatively simple and 
20 easily optimized and the transformation process is 

relatively fast. If the device profile link 412 is a look- 
up table, the device link calculator 414 optionally uses 

20 



linear interpolation to refine the destination device 



coordinates . 



The device link calculator 414 can be 



implemented, e . gr 



• / 



using a conventional multidimensional 



linear interpolator . 
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FIG. 5 illustrates an example device link 



generator 500 that implements part of the color transformer 
400 of FIG. 4. The device link generator 500 includes a 
device link table builder 502 that creates a look-up table 
to enable rapid interpolation of destination device 
10 coordinates from source device coordinates. It should be 

understood that if the device profile link is a mathematical 
expression rather than as a look-up table, an analogous 
transformation generator replaces the device link table 
builder 502. Such a transformation builder may, for 
- 15 example, generate coefficients for use in the mathematical 
expression. To facilitate the discussion, however, the 
device link generator 500 is assumed to include a device 
link table builder 502. The device link table builder 502 
generates the look-up table by generating a series of source 
20 device coordinates as input value entries and determining 

the optimal destination device coordinates as output values 
corresponding to the input values. The device link table 
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builder 502 generates all combinations of source device 
coordinates using, for example, a series of nested loops, 
one loop for each dimension of the source device color 
coordinate space . 
5 To reduce the computational and memory 

requirements for constructing and storing the look-up table, 
the look-up table typically contains a relatively small 
number of entries along each dimension. With a relatively 
small table, interpolation is used to convert source 

10 coordinates to destination coordinates. The total number of 
entries in the look-up table can be expressed as D d N s d , where 
d is the dimensionality of the source device color space, D d 
is the dimensionality of the destination device color space, 
and N s is the number of entries along each dimension of the 

15 look-up table. For example, a look-up table that is used to 
transform color coordinates between two CMYK (i.e., four- 
dimensional) color spaces can contain 4xl7 4 , or 334,084 
entries . 

It should be understood that the look-up table 
2 0 need not have the same number of entries along each 

dimension. If the look-up table contains N k entries along 
each respective dimension, where k ranges from 1 to d, the 
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total number of entries in the look-up table can be 

d 

expressed as i^FI^A - For example, a look-up table that is 

used to transform color coordinates between two CMYK color 
spaces and that is to have fifteen entries along three 
5 dimensions and seventeen entries along one dimension 
contains 4x15x15x15x17, or 229,500 entries. 

In application environments in which it is 
desirable to further reduce computational and memory 
requirements, the device link table builder 502 may select 
10 only a subset of the total number of entries in each 

dimension of the look-up table, perform the method loop 
calculations using that subset, and perform, for example, a 
spline interpolation to fill in the remaining entries of the 
look-up table. 

15 The device link table builder 502 provides PCS 

parameters and source device coordinates to a source device 
profile interpreter 504. The source device profile 
interpreter 504 generates source PCS values and provides the 
source PCS values and the source device coordinates to an 

2 0 error reducer 506. In a specific embodiment, the error 
reducer 506 is implemented using an error minimization 
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technique. Using the source device coordinates, the error 
reducer 506 estimates an initial set of destination 
coordinates that are likely to realize an accurate color 
match with the source device coordinates. This estimation 
5 process may be performed using a relatively simple 
technique. For example, for estimating destination 
coordinates in an RGB space corresponding to source 
coordinates in a CMYK color space, the estimation process 
may use the following equations: 
10 C = 1 - R 

M = 1 - G 

Y = 1 - B 

Alternatively, the source coordinates can be used to 
estimate the destination coordinates if the source and 
15 destination imaging systems use similar color coordinate 
spaces . 

The error reducer 506 provides the set of 
estimated destination device coordinates to a destination 
device profile interpreter 508, which also receives the PCS 
20 parameters from the device link table builder 502. The 

destination device profile interpreter 508 then generates a 
set of destination PCS values as a function of the estimated 
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destination device coordinates and the PCS parameters and 
provides the destination PCS values to the error reducer 
506. If the error between the destination PCS values and 
the source PCS values is non-zero, the error reducer 506 
5 uses an error reduction (e.g., an error minimization) 
technique to reduce the error between the source and 
destination PCS values. In one embodiment, this is 
implemented by repeatedly querying the destination device 
profile interpreter 508 with selected estimates of 

10 destination device coordinates. This process can continue 
until destination device coordinates are found that satisfy 
a quality threshold, for example, that yield the minimum 
error. The error reducer 506 returns these destination 
device coordinates to the device link table builder 502, 

15 which enters them in an appropriate location in the look-up 
table. The device link table builder 502 then enters the 
next set of table input entries corresponding to a set of 
source device color coordinates . 

For colors within the gamut of the destination 

20 device, the error can be reduced using any of a variety of 
reduction techniques. For example, Powell's method can be 
used to perform the error reduction or error minimization. 
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For additional information regarding Powell's method, 
reference may be made to William H. Press et al . , Numerical 
Recipes in C (1992), pp. 309-315, available from Cambridge 
University Press, incorporated herein by reference and 
5 attached hereto as Appendix A. 

Using this approach, the error reducer 506 
generally defines an error function having input parameters 
that can be varied by the error reduction technique. The 
error reducer 506 then determines the optimal values of the 
10 input parameters resulting in a minimal error. To determine 
the values of destination device coordinates using the 
minimum PCS error between the source and destination PCS 
values, the variable input parameters are the destination 
device coordinates. Accordingly, in this specific 
> 15 implementation, the error reducer 506 defines the error 
function as: 

Error (D) = AE (R s , R d (D) ) 
where D is a vector defined by the destination device 
coordinates, R s is a vector defined by the source PCS values 
20 and R d (D) is a vector function producing destination PCS 

values as a function of the destination coordinate vector D, 
and AE is the Euclidean distance error between R s and R d (D) . 
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The Euclidean distance error may be expressed using the 
following equation: 

AE(R 1# R 2 ) = ^(L, - ij 2 + (a, - a 2 ) 2 + (h L ~ b 2 f 
The above equation assumes that the PCS is 
5 implemented as the CIELAB color space. It should be 

understood, however, that other color spaces may be used as 
a PCS. For example, one color space that is particularly 
suited for use as a PCS is described in the previously- 
referenced copending U.S. Patent Application, entitled 
10 "Characterization of Color Imaging Systems." 

Using this same approach, a non-zero optimal error 
indicates that the source device is out of gamut relative to 
the destination device at that location in the PCS. In such 
situations, the error reducer 506 optionally uses the 
15 destination device coordinates that result in a minimum Ae 
value. Alternatively, the error reducer 506 may use these 
values as an initial estimate and recalculate the optimal 
destination device coordinates using a new error function 
that employs weighting factors, polar coordinates in the 
2 0 chromatic plane of the PCS space, or both. 

The error reducer 506 optionally uses a gamut 
mapping parameter received from the device link table 
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builder 502 to decide how to map coordinates that are out of 
gamut relative to the source device. For example, the gamut 
mapping parameter may specify modes in which each technique 
is used for obtaining destination device coordinates. One 
mode, for example, may use lightness, chroma, and hue values 
L , C , and h instead of LAB values L , a , and b , where: 




Another mode uses the above lightness, chroma, and hue 
values as well as weighting factors: 

Error(D) = AEW(R S , R d (D) , W) 

Aew(r 1 , r 2/ w) = a /w l (l l - l 2 ) 2 + w c (q - c 2 ) 2 + - h^f 
where the PCS vectors R s , R d (D) are converted to lightness, 
chroma, and hue values either before or after passing them 
to the error function. 

If the weighting factors are one, the above- 
weighted error reduction function AEW() is identical to the 
standard AE() error reduction function. It should be 
noted, however, that weighting factors of W L =3 , W c =l, and 
W h =1.5 yield particularly accurate visual results. These 
weighting factors produce an error function that gives 
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priority first to lightness, then to hue, then to chroma. 
These weighting factors can also be provided to the error 
reducer 506 as gamut mapping parameters by the device link 
table builder 502. 

Creating the device profile link via error 
reduction of the forward transformations of the devices 
realizes a number of advantages. For example, errors in 
color conversion are limited to those attributable to 
rounding and interpolation. As a result, the cumulative 
error from repeated forward and reverse transformations 
between the source and destination device color spaces is 
substantially reduced. Additionally, the color transformer 
can select the gamut mapping technique . The color 
transformer can rely on the forward transform information 
and realize consistent gamut mapping between device profiles 
supplied by different vendors. It should be noted that 
errors due to interpolation of the device profile link 
decrease as the number of table entries in each dimension of 
the look-up table increases toward the maximum number of 
gray levels. This error also decreases if a one-dimensional 
tone reproduction table is used to transform the color 
values. For additional information regarding the use of a 



one-dimensional tone reproduction table, reference is made 
to U.S. Patent No. 5,432,906, issued to Gary H. Newman, 
assigned to Eastman Kodak Company, and incorporated by 
reference . 

5 Creating the device profile link using error 

reduction also allows transformation between CMYK device 
spaces that maps the tone response of the source and 
destination black (K) channels while maintaining an accurate 
match with the L*a*b* data. For transformation from an RGB 

10 source device to a CMYK destination device, the RGB color 
coordinates used by the source device lack K channel 
information . Some conventional color transformation 
techniques use a process known as gray component removal 
(GCR) to define a relationship between K values and CMY 

15 values in the reverse transformation (i.e., L*a*b* to CMYK). 
For example, the reverse transformation may be performed 
with K initially set to zero. The value of K can then be 
calculated based on the minimum of the C, M, and Y values. 
The CMY values can then be recalculated using an algebraic 

20 calculation or using the forward model to obtain the closest 
value of L*a*b* input using the new calculated K value. This 
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process involves a reverse transformation from Lab color 
values to CMYK color values with a fixed definition of GCR. 

This process, however, loses the K channel 
information or the CMY channel information during the 
translation between CMYK color spaces because the source 
color values are transformed to a three-dimensional 
intermediary color space during conversion to destination 
CMYK values. To preserve the K channel information, the 
error reducer 506 determines optimal K values in the 
destination color space that correspond to the K values in 
the source device color space, e.g., values between 0 and 
2 55. These values can be created, for example, by 
generating a series of source K values ranging from minimum 
to maximum, fixing the source and destination CMY values at 
0, and finding destination K values with minimum AE error 
relative to each of the source K values. These source and 
destination K values can be loaded into a lookup table for 
quick conversion of source K to destination K values. By 
using error reduction to determine optimal K values in the 
destination color space, the device link generator 500 
preserves K channel information. This results in improved 



accuracy of the K channel information when converting colors 
between CMYK devices. 

After loading the source and destination K values 
into a lookup table, when the error reducer 506 receives 
5 source L*a*b* and CMYK values, the error reducer 506 

initially maps the source K channel to the destination K 
channel. The error reduction procedure is then used for 
varying the destination CMY values to obtain the best match 
for the respective L*a*b* values. If Ae = 0, control 

10 returns to the device link table builder 502, which enters 
the calculated destination CMYK values into the device link 
table. If Ae is greater than zero, then the destination 
CMY values corresponding to the destination K value in 
question are out of gamut relative to the target L*a*b* 

15 values. This may be, for example, because the source CMY 

values corresponding to K = 0 result in a color that is out 
of gamut with the destination device, or because the 
destination K value in the particular region of destination 
CMY color space is either too high or too low, i.e., the 

2 0 mixture of K with CMY is such that the resulting color is 

too dark or too light relative to the targeted L*a*b* value. 



32 



To reduce the AE error, K can be varied in a 

* * * 

controlled way so as to ensure both optimal Lab color and 
optimal matching of the K source channel behavior. This can 
be performed, for example, by alternately fixing the current 
CMY values while performing error reduction on variable K 
values and fixing the K value while performing error 
reduction on variable CMY values. When it is determined 
that neither varying CMY nor varying K improves the AE 
error, it can be assumed that the optimal CMYK values have 
been determined to satisfy both the color matching and K 
channel accuracy criteria. Control then optionally returns 
to the device link table builder 502. While the above 
discussion assumes that the error reducer 506 performs the 
mapping between source and destination K values, it should 
be understood that the device link table builder 502 can 
perform the mapping. 

It should be understood that other approaches can 
be used to improve the accuracy of the K channel 
information. For example, the PCS can be implemented as a 
color space having the same number of dimensions, e.g., 
four, as CMYK space. Using a PCS having the same 
dimensionality as the device space prevents the loss of 



color channel information. In a specific example 
embodiment, the first three channels of this PCS are the PCS 
currently used by the system {e.g., LAB, L*a*b*, or XYZ) . 
The fourth channel indicates a PCS value indicative of the 
black channel or relating to the black channel (e.gr., L* or 
tristimulus value Y) . The process can be performed in a 
manner similar to that performed by the ICC specification as 
in, for example, ColorSync 2.1 available from Apple 
Computer. 

FIG. 6 illustrates an example color transformation 
method 60 0 according to the present invention. At a block 
602, selected source device color coordinates are mapped to 
a PCS. Destination device color coordinates are then 
estimated as a function of the source device color 
coordinates, as depicted at a block 6 04. These estimated 
destination device color coordinates are then mapped to the 
PCS at a block 606 . 

At a block 6 08, an error between the PCS values 
corresponding to the source and destination device color 
coordinates is determined. At a decision block 610, the 
method determines whether the error satisfies a quality 
criterion, such as error minimization. In certain 



applications, the quality criterion can be defined as 
reduction of the error below a threshold value. If the 
error does not satisfy the quality criterion, flow proceeds 
to a block 612, at which the estimated destination device 
5 color coordinates are adjusted to reduce the error. This 
process repeats until the error is reduced. 

After the error is reduced, flow proceeds to a 
block 614, at which the optimal destination device color 
coordinates thus obtained are entered into a color map. 

10 Next, the method determines whether the color map is filled, 
as depicted at a decision block 616. If the color map 
contains empty entries, flow proceeds to a block 618. New 
source device color coordinates are then selected, and then 
flow returns to the block 602. This process continues until 

15 the color map is filled. The color map can then be stored 

as, for example, a data file for future reference. The user 
can specify the desired source, destination, and 
intermediate profiles and the user preferences used to 
generate the device profile link. Upon recognizing that a 

20 color map has already been developed for a particular 

combination of these profiles, the system can load the data 
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file. Loading the data file instead of reconstructing the 
color map saves computation time and other resources. 

The device profile link can be generated each time 
the user requests a new combination of device profiles. 
5 Alternatively, the user can specify in advance a series of 
source, intermediate, and destination profiles and allow the 
system to preprocess these lists of profiles into their 
respective device profile links and store them. When the 
user requests that a particular transform be performed on 

10 image data using a previously defined combination of source, 
intermediate, and destination profiles, the system retrieves 
the associated device profile link. Retrieving the device 
profile link improves the processing speed. 

While the above discussion has assumed that the 

15 device profile link describes a conversion between two 

device profiles, it should be understood that the device 
profile link can be used to describe a conversion between 
any number of device profiles. For example, N device 
profiles can be concatenated using a single device profile 

2 0 link. To concatenate the device profiles, the color 

conversion is performed using the PCS to convert colors 
between each device profile to be concatenated. Performing 
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error reduction on the forward transforms between the 
individual device profiles improves the accuracy of the 
concatenated device profile link between the first and n th 
device profiles. 
5 The various embodiments described above are 

provided by way of illustration only and should not be 
construed to limit the invention. Those skilled in the art 
will readily recognize various modifications and changes 
that may be made to the present invention without strictly 
10 following the example embodiments and applications 

illustrated and described herein, and without departing from 
the true spirit and scope of the present invention, which is 
set forth in the following claims. 



15 What is claimed is; 
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1 l. For use in transforming colors between color 

2 imaging systems, a color mapping method comprising: 

3 using forward transformation profiles that 

4 characterize the color imaging systems to generate 

5 respective sets of device -independent color values for the 

6 color imaging systems; 

7 calculating color conversions by recursively 

8 reducing differences between the sets of device -independent 

9 color values; and 

10 constructing a color map describing a relationship 

11 between the color imaging systems using the color 

12 conversions . 

1 2. A color mapping method, according to claim 1, 

2 further comprising recursively reducing differences between 

3 black channel information. 

1 3. A color mapping method, according to claim 1, 

2 further comprising using an error function for calculating 

3 the color conversions. 

1 4. A color mapping method, according to claim 1, 

2 further comprising configuring at least one of the profiles 
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3 to account for certain perceptual effects on color 

4 appearance . 

1 5. A color mapping method, according to claim 1, 

2 wherein the color map comprises at least one of the 

3 following: a lookup table 7 and an equation. 

1 6. A color mapping method, according to claim 1, 

2 further comprising: 

3 storing the color map; 

4 detecting respective types of color imaging 

5 devices between which a color transformation is to be 

6 performed; and 

7 in response to the detected types, selecting a 

8 stored color map. 

1 7 . For use in transforming colors between source 

2 and destination color imaging systems, a color mapping 

3 method comprising : 

4 using profiles that characterize the color imaging 

5 systems to generate device- independent color values for the 

6 source color imaging system, the device -independent color 
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7 values having a same dimensionality as the source color 

8 imaging system; 

9 using the profiles to perform a color conversion 

10 for converting the device -independent color values to 

11 device -dependent values of the destination color imaging 

12 system; and 

13 using the color conversion to define a color map 

14 for transforming colors between the color imaging systems. 

1 8. A color mapping method, according to claim 7, 

2 wherein the color conversion is performed at least twice. 

1 9. A color mapping method, according to claim 7, 

2 further comprising: 

3 using the color conversion to evaluate its 

4 accuracy at least once; and 

5 revising the color conversion at least once to 

6 improve its accuracy. 

1 10. For use in transforming colors between source 

2 and destination color imaging systems, a color mapping 

3 method comprising: 
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1 (a) using profiles characterizing the color 

2 imaging systems to generate device -independent color values 

3 for the source color imaging system, the device -independent 

4 color values having a same dimensionality as the source 

5 color imaging system; 

6 (b) using the profiles to perform a color 

7 conversion for converting the device- independent color 

8 values to device -dependent values of the destination color 

9 imaging system; 

10 (c) using the color conversion to improve the 

11 accuracy of the color conversion relative to a quality 

12 threshold; 

13 (d) returning to step (c) until the color 

14 conversion satisfies the quality threshold; and 

15 (e) using the color conversion to define a color 

16 map for transforming colors between the color imaging 

17 systems . 

1 11. For use in transforming colors between color 

2 imaging systems, a color mapping arrangement comprising: 

3 means for using forward transformation profiles 

4 that characterize the color imaging systems to generate 



41 



5 respective sets of device- independent color values for the 

6 color imaging systems; 

7 means for calculating color conversions by 

8 recursively reducing differences between the sets of device- 

9 independent color values; and 

10 means for constructing a color map describing a 

11 relationship between the color imaging systems using the 

12 color conversions . 

1 12. For use in transforming colors between first 

2 and second color imaging systems respectively using first 

3 and second color coordinate systems, a color mapping method 

4 comprising : 

5 (a) generating first device- independent color 

6 coordinates as a function of color coordinates in the first 

7 color coordinate system; 

8 (b) estimating preliminary color coordinates in 

9 the second color coordinate system; 

10 (c) generating second device- independent color 

11 coordinates as a function of the preliminary color 

12 coordinates ; 
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13 (d) adjusting the preliminary color coordinates 

14 to reduce an error between the first and second device- 

15 independent color coordinates; 

16 (e) returning to step (a) until the error 

17 satisfies a quality threshold; and 

18 (f) constructing a color map describing a 

19 relationship between the first and second color imaging 

20 systems as a function of the adjusted color coordinates. 

1 13 . A color mapping method, according to claim 

2 12 , further comprising using the color coordinates in the 

3 first color coordinate system to estimate the preliminary 

4 color coordinates. 

1 14 . For use in transforming colors between color 

2 imaging systems , a color mapping arrangement comprising: 

3 a computer arrangement, programmed to 

4 use forward transformation profiles that 

5 characterize the color imaging systems to generate 

6 respective sets of device -independent color values for the 

7 color imaging systems, 
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8 calculate color conversions by recursively 

9 reducing differences between the sets of device -independent 

10 color values, and 

11 construct a color map describing a 

12 relationship between the color imaging systems using the 

13 color conversions; and 

14 a memory, configured and arranged to store the 

15 color map. 

1 15. A color mapping arrangement, according to 

2 claim 14, wherein the computer arrangement is further 

3 programmed to use an error function for calculating the 

4 color conversions. 

1 16. A color mapping arrangement, according to 

2 claim 14, wherein the computer arrangement is further 

3 programmed to configure at least one of the profiles to 

4 account for certain perceptual effects on color appearance. 

1 17. A color mapping arrangement, according to 

2 claim 14, wherein the computer arrangement is further 

3 programmed to construct at least one of the following: a 

4 lookup table, and an equation. 
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1 18. A color mapping arrangement, according to 

2 claim 14, wherein the computer arrangement is further 

3 programmed to 

4 detect respective types of color imaging devices 

5 between which a color transformation is to be performed, and 

6 in response to the detected types, select a stored 

7 color map . 

1 19. For use in transforming colors between color 

2 imaging systems, a data storage medium storing a computer- 

3 executable program that, when executed, 

4 uses forward transformation profiles that 

5 characterize the color imaging systems to generate 

6 respective sets of device -independent color values for the 

7 color imaging systems; 

8 calculates color conversions by recursively 

9 reducing differences between the sets of device- independent 

10 color values, and 

11 constructs a color map describing a relationship 

12 between the color imaging systems using the color 

13 conversions . 
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1 20, A data storage medium, according to claim 19, 

2 wherein the computer-executable program recursively reduces 

3 differences between black channel information. 

1 21. A data storage medium, according to claim 19, 

2 wherein the computer- executable program uses an error 

3 function for calculating the color conversions. 

1 22. A data storage medium, according to claim 19, 

2 wherein the computer-executable program configures at least 

3 one of the profiles to account for certain perceptual 

4 effects on color appearance. 



1 23 . A data storage medium, according to claim 19, 

2 wherein the computer-executable program generates at least 

3 one of the following: a lookup table, and an equation. 

1 24 . A data storage medium, according to claim 19, 

2 wherein the computer-executable program: 

3 stores the color map; 

4 detects respective types of color imaging devices 

5 between which a color transformation is to be performed; and 
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6 in response to the detected types, selects a 

7 stored color map. 
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Abstract 



A color mapping method is used in transforming 



colors between color imaging systems . The method includes 
using forward transformation profiles that characterize the 
5 color imaging systems to generate respective sets of device- 
independent color values for the color imaging systems. 
Color conversions are calculated by recursively reducing 
differences between the respective sets of device- 
independent color values. Based on these color conversions, 
10 a color map is constructed that describes a relationship 
between the color imaging systems . 
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We have just seen that choosing N points uniformly randomly in an n- b€?^o' 

dimensional space leads to an error term in Monte Carlo integration that decreases § IF 

as l/VN, In essence, each new point sampled adds linearly to an accumulated sum o JL 1 1 1- 

that will become the function average, and also linearly to an accumulated sum of g s* = K 3 
squares that will become the variance (equation 7.6.2). The estimated error comes 3 *© c 

from the square root of this variance, hence the power N* */ 2 . f 1 3 i S 

Just because this square root convergence is familiar does not, however, mean j§ jj % | g 

that it is inevitable. A simple counterexample is to choose sample points that lie i s § % Z 

on a Cartesian grid, and to sample each grid point exactly once (in whatever order). " f 2 § Q 

The Monte Carlo method thus becomes a deterministic quadrature scheme — albeit | -Z. * ^, S 

a simple one — whose fractional error decreases at least as fast as N * 1 (even faster > §--! <| s 

• 3 3^ S r> 

if the function goes to zero smoothly at the boundaries of the sampled region, or © *< s* | - 

is periodic in the region). | Hie? In 

The trouble with a grid is that one has to decide in advance how fine it should .§; §o 1.5 

be. One is then committed to completing all of its sample points. With a grid, it is o £ * 2 o 

not convenient to "sample untiF some convergence or termination criterion is met. 8 <> | -2 « 

One might ask if there is not some intermediate scheme, some way to pick sample * | » g g 

points "at random," yet spread out in some self-avoiding way, avoiding the chance 1 1 1 <S ^ 

clustering that occurs with uniformly random points. ? I ? ^ o 

A similar question arises for tasks other than Monte Carlo integration. We might | S. | ^ § 

want to search an n-dimensional space for a point where some (locally computable) ® g 3 f c 

condition holds. Of course, for the task to be computationally meaningful, there v "g <| g =? 

had better be continuity, so that the desired condition will hold in some finite n- I g S 

dimensional neighborhood. We may not know a priori how large that neighborhood - " 9 — 



is, however. We want to "sample untir the desired point is found, moving smoothly * 1 2 % m 
to finer scales with increasing samples. Is there any way to do this that is better 1 1-< |S 
than uncorrelated, random samples? zol P s 

The answer to the above question is "yes." Sequences of n-tuples that fill §. o<o g 
n-space more uniformly than uncorrelated random points are called quasi-random > § f ~ 
sequences. That term is somewhat of a misnomer, since there is nothing "random" ® | 
about quasi-random sequences: They are cleverly crafted to be, in fact, swfc-random. II 3 
The sample points in a quasi-random sequence are, in a precise sense, "maximally 
avoiding" of each other. 

A conceptually simple example is Halton's sequence [1 ]. In one dimension, the 
jth number Hj in the sequence is obtained by the following steps: (i) Write j as a 
number in base &, where b is some prime. (For example j = 17 in base b = 3 is 
122.) (ii) Reverse the digits and put a radix point (i.e., a decimal point base 6) in 
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points i to 128 




.2 .4 .6 .8 
points 513 to 1024 




0 .2 .4 .6 .8 1 
points 129 to 512 




.2 .4 .6 .8 
points 1 to 1024 



figure 7.7.1. first 1024 points of a two-dimensional Sobol' sequence. The sequence is generated 
number-theoretically* rather than randomly, so successive points at any stage "know" how to fill in the 
gaps in the previously generated distribution. 



front of the sequence. (In the example, we get 0.221 base 3.) The result is Hj. To 
get a sequence of n-tuples in n-space, you make each component a Halton sequence 
with a different prime base b. Typically, the first n primes are used. 

It is not hard to see how Halton's sequence works: Every time the number of 
digits in j increases by one place, j's digit-reversed fraction becomes a factor of 
b finer-meshed. Thus the process is one of filling in all the points on a sequence 
of finer and finer Cartesian grids — and in a kind of maximally spread-out order 
on each grid (since, e.g., the most rapidly changing digit in j controls the most 
significant digit of the fraction). 

Other ways of generating quasi-random sequences have been suggested by 
Faure, Sobol*, Niederreiter, and others. Bratley and Fox [2] provide a good review 
and references, and discuss a particularly efficient variant of the Sobol' [3] sequence 
suggested by Antonov and Saleev [4]. It is this Antonov-Saleev variant whose 
implementation we now discuss. 
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Degree 


Primitive Polynomials Modulo 2* 


I 


0 (i.e., x + 1) 


2 


I (i.e., x 2 + x + 1) 


3 


1, 2 (i.e„ x 3 + x + 1 and x 3 4- x 2 + 1) 


4 


1, 4 (i.e„ x 4 4* x + 1 and x 4 + x 3 + 1) 


c 
J 




6 


1, 13, 16, 19, 22, 25 


7 


1, 4, 7. 8, 14, 19, 21, 28, 31, 32, 37, 41, 42, 50, 55. 56, 59, 62 


8 


14, 21, 22, 38, 47, 49, 50, 52, 56, 67, 70, 84, 97, 103, 1 15, 122 


9 


8, 13, 16, 22, 25, 44, 47, 52, 55, 59, 62, 67, 74, 81, 82, 87, 91, 94, 103, 104, 109, 122, 
124, 137, 138, 143, 145, 152, 157, 167, 173, 176, 181, 182, 185, 191, 194, 199,218,220, 
227, 229, 230, 234, 236. 241, 244, 253 


10 


4, 13, 19, 22, 50, 55, 64, 69, 98, 107, 115, 121, 127, 134, 140, 145, 152, 158, 161, 171, 
181, 194, 199, 203, 208, 227, 242, 251, 253, 265, 266, 274, 283, 289, 295, 301, 316, 
319, 324, 346, 352, 361, 367, 382, 395, 398, 400, 412, 419, 422, 426, 428, 433, 446, 
454, 457, 472, 493, 505, 508 


♦Expressed as a decimal integer representing the interior bits (that is, omitting the 
high-order bit and the unit bit). 



The Sobol' sequence generates numbers between zero and one directly as binary fractions 
of length w bits, from a set of w special binary fractions, V5, i = 1, 2, — ,ti;, called direction 
numbers. In SoboPs original method, the jth number Xj is generated by XORing (bitwise 
exclusive or) together the set of Vi 's satisfying the criterion on i, "the ith bit of j is nonzero." 
As j increments, in other words, different ones of the Vi*s flash in and out of Xj on different 
time scales. V\ alternates between being present and absent most quickly, while Vk goes from 
present to absent (or vice versa) only every 2 fc 1 steps. 

Antonov and Saleev's contribution was to show that instead of using the bits of the 
integer j to select direction numbers, one could just as well use the bits of the Gray code of j 9 
G(j). (For a quick review of Gray codes, look at §20.2.) 

Now G(j) and G(J 4- 1) differ in exactly one bit position, namely in the position of the 
rightmost zero bit in the binary representation of j (adding a leading zero to j if necessary). A 
consequence is that the j + 1st SoboF-Antonov-Saleev number can be obtained from the jth 
by XORing it with a single Vi, namely with i the position of the rightmost zero bit in j. This 
makes the calculation of the sequence very efficient, as we shall see. 

Figure 7.7. 1 plots the first 1024 points generated by a two-dimensional Sobol' sequence. 
One sees that successive points do "know" about the gaps left previously, and keep filling 
them in, hierarchically. 

We have deferred to this point a discussion of how the direction numbers Vi are generated. 
Some nontrivial mathematics is involved in that, so we will content ourself with a cookbook 
summary only: Each different Sobol' sequence (or component of an n-dimensional sequence) 
is based on a different primitive polynomial over the integers modulo 2, that is, a polynomial 
whose coefficients are either 0 or 1, and which cannot be factored (using modulo 2 integer 
arithmetic) into polynomials of lower order. (Primitive polynomials modulo 2 were used in 
§7.4, and are further discussed in §20.3.) Suppose P is such a polynomial, of degree 
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Initializing Values Used in sobseq 


Degree 


Polynomial 


Starting Values 


I 


0 




(3) 


(5) 


05) ... 


2 


1 




1 


(7) 


(11) ... 


3 


1 




3 


7 


(5) ... 


3 


2 




3 


3 


05) ... 


4 


1 




1 


3 


13 ... 


4 


4 




1 


5 


9 ... 


Parenthesized values are not freely specifiable, but are forced by the required recurrence 
for this degree. 



Define a sequence of integers Mi by the ?-term recurrence relation, 

Mi = 2ai Mi. i 0 2 2 a 2 Mi. 2 0 • • • 0 2* x Mi. q + x a q . 1 0 (2*JWi. , 0 M». q ) (7.7.2) 

Here bitwise XOR is denoted by 0. The starting values for this recurrence are that Mi , . . . , M q 
can be arbitrary odd integers less than 2, . . . , 2 9 , respectively. Then, the direction numbers 
Vi are given by 



V i = M i /2 i i=l,. 



(7.7.3) 



The accompanying table lists all primitive polynomials modulo 2 with degree q < 10. 
Since the coefficients are either 0 or 1 , and since the coefficients of x q and of 1 are predictably 
1, it is convenient to denote a polynomial by its middle coefficients taken as the bits of a binary 
number (higher powers of x being more significant bits). The table uses this convention. 

Turn now to the implementation of the Sobol' sequence. Successive calls to the function 
sobseq (after a preliminary initializing call) return successive points in an n-dimensional 
Sobor sequence based on the first n primitive polynomials in the table. As given, the 
routine is initialized for maximum n of 6 dimensions, and for a word length w of 30 bits. 
These parameters can be altered by changing MAXBIT (= w) and MAXDIM, and by adding 
more initializing data to the arrays ip (the primitive polynomials from the table), mdeg (their 
degrees), and iv (the starting values for the recurrence, equation 7.7.2). A second table, 
above, elucidates the initializing data in the routine. 

#include "nrutil.h** 
#def ine MAXBIT 30 
#def ine HAXDIH 6 

void sobseqCint *n, float x[3) 

When n. is negative, internally initializes a set of MAXBIT direction numbers for each of MAXDIM 
different Sobol r sequences. When n is positive (but <MAXDIM), returns as the vector x[l . .n] 
the next values from n of these sequences, (n must not be changed between initializations.) 
{ 

int j,k,l; 

unsigned long i t isi,ipp; 
static float fac; 

static unsigned long in,ix[MAXDIM+l] ,*iu[MAXBIT+i] ; 
static unsigned long mdeg[MAXDIM+i]*<0,l,2,3,3,4,4}; 
static unsigned long ip[MAXDIM+i3-{0,0,i,i,2,l,4}; 
static unsigned long iv [MAXDIM*MAXBIT+l3 -{ 

0,1,1,1,1,1,1,3,1,3,3,1,1,5,7,7,3,3,5,15,11,5,15,13,9}; 
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if (*n < 0) { 

for (k»l;k<»KAXDIM;k++) ix[k>0; 



Initialize, don't return a vector. 
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in-O; 

if (ivCi] !- 1) rttura; 
fac-1.0/(lL « MAXBIT) ; 

lor <j-l,k-0;j<-MAXBIT;j++,k+-MAXDIM) iu[j] « *iv[k] ; 
To allow both ID and 2D addressing, 
for (k-l;k<«MAXDIM;k++) { 

for <j-l;j<«*d«gCk];j++> iuCj] Ck] «- (MAXBIT- j) ; 

Stored values only require normalization. 

for < j-md«g[k]+l; j<-MAXBIT; j++) { Use the recurrence to get other val- 
ipp-ipCk] ; ues. 
i-iutj-mdegDc]][k]; 
i (i » mdegCk]); 
for (l-Bdeg[k]-l;l>«i;l--) { 

if (ipp ft 1) i ~- iuCj-ilCk]; 

ipp »- i; 

> 

iuEj3 Ck]-i; 



///« 



ST o : 

-r CO ' 



o 3. < 



y 

y else { Calculate the next vector in the se- 

im-in++; quence. 

for (j«i; j<-HAXBIT; j++) { Find the rightmost zero bit. £ j 

if (Kim fcl)) break; 
im »• i; 

if (j > MAXBIT) nrerrorCMAXBIT too small in sobseq"); fe-§ 3 3»5 

im-(j-i)*MAXDIM; "IS i Q 

for (k«i;k<»IMIN(*n,HAXDIM);k++) -C XOR the appropriate direction num- o~2% 

ix[k] "« ivEim+k] ; ber into each component of the §r £ 8 3 52 

x [k] ~ix Ck] *f ac ; vector and convert to a floating 

} number. _ 

> 8 s SQ] 



How good is a Sobol' sequence, anyway? For Monte Carlo integration of a smooth 8 ° 1 § UJ 

function in n dimensions, the answer is that the fractional error will decrease with N, the = |[ o ^ o 

number of samples, as (In N) n /N t i.e., almost as fast as 1 /N. As an example, let us integrate | 5 £. g g 

a function that is nonzero inside a torus (doughnut) in three-dimensional space. If the major - = I g §j 

radius of the torus is Ho, the minor radial coordinate r is defined by ° I -n o 

r = ([(x 2 + y*?» - ilo] 2 + * 2 ) V2 (7.7.4) 1 1 1 1 | 

Let us try the function S f. c g 0 

H X ,V,z) = \l+«*(-s) r< n (7 . 7 . 5) | t jl 

[0 r>r 0 sS-llt 

which can be integrated analytically in cylindrical coordinates, giving z o •< P g 

o 7" ~ 1 

dx dy dz f{x, y, z) = 2n 2 a 2 Ro (7.7.6) > \ 



With parameters Rq — 0.6, ro = 0.3, we did 100 successive Monte Carlo integrations of A ? 
equation (7.7.4), sampling uniformly in the region -1 < x,y,z < 1, for the two cases of 
uncorrected random points and the Sobol' sequence generated by the routine sobseq. Figure 
7.7.2 shows the results, plotting the r.m.s. average error of the 100 integrations as a function 
of the number of points sampled. (For any single integration, the error of course wanders 
from positive to negative, or vice versa, so a logarithmic plot of fractional error is not very 
informative.) The thin, dashed curve corresponds to uncorrected random points and shows 
the familiar iV 1/2 asymptotics. The thin, solid gray curve shows the result for the SoboF 



( 



314 



Chapter 7. Random Numbers 



B 
.5 



u 

•a 

§ 

•■a 

1 



.01 r 



.001 



i i i i I 



"i 1 — i — i — i — rn 



-T 1 — i — i — i — PT 



•ill * % 






pseudo-random, hard boundary ■ ■■■ 

pseudo-random, soft boundary 

quasi-random, hard boundary ^™ 
quasi-random, soft boundary 



-J 1 I 1 1 — L 



100 



1000 10000 
number of points N 



105 



Figure 7.7.2. fractional accuracy of Monte Carlo integrations as a function of number of points sampled, 
for two different integrands and two different methods of choosing random points. The quasi-random 
Sobol' sequence converges much more rapidly than a conventional pseudo-random sequence. Quasi- 
random sampling does better when the integrand is smooth ("soft boundary") than when it has step 
discontinuities ("hard boundary"). The curves shown are the r.m.s. average of 100 trials. 

sequence. The logarithmic term in the expected (In N) 3 /N is readily apparent as curvature 
in the curve, but the asymptotic N' 1 is unmistakable. 

To understand the importance of Figure 7.7.2, suppose that a Monte Carlo integration of 
/ with 1% accuracy is desired. The Sobol' sequence achieves this accuracy in a few thousand 
samples, while pseudorandom sampling requires nearly 100,000 samples. The ratio would 
be even greater for higher desired accuracies. 

A different, not quite so favorable, case occurs when the function being integrated has 
hard (discontinuous) boundaries inside the sampling region, for example the function that is 
one inside the torus, zero outside, 



/(z>2/>*) = {o 



(7.7.7) 



r < ro 
r > ro 

where r is defined in equation (7.7.4). Not by coincidence, this function has the same analytic 
integral as the function of equation (7.7.5), namely 2n 2 a 2 Ro. 

The carefully hierarchical Sobol' sequence is based on a set of Cartesian grids, but the 
boundary of the torus has no particular relation to those grids. The result is that it is essentially 
random whether sampled points in a thin layer at the surface of the torus, containing on the 
order of N 2 ^ 3 points, come out to be inside, or outside, the torus. The square root law, applied 
to this thin layer, gives iV 1/3 fluctuations in the sum, or N' 2/3 fractional error in the Monte 
Carlo integral. One sees this behavior verified in Figure 7.7.2 by the thicker gray curve. The 
thicker dashed curve in Figure 7.7.2 is the result of integrating the function of equation (7.7.7) 
using independent random points. While the advantage of the Sobo!' sequence is not quite so 
dramatic as in the case of a smooth function, it can nonetheless be a significant factor (~5) 
even at modest accuracies like 1%, and greater at higher accuracies. 
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Note that we have not provided the routine aobaeq with a means of starting the 
sequence at a point other than the beginning, but this feature would be easy to add. Once 
the initialization of the direction numbers iv has been done, the jth point can be obtained 
directly by XORing together those direction numbers corresponding to nonzero bits in the 
Gray code of j, as described above. 

The Latin Hypercube S3 t,os 

We might here give passing mention the unrelated technique of Latin square or I » § ^ i 

Latin hypercube sampling, which is useful when you must sample an iV-dimensional * £ | 

space exceedingly sparsely, at M points. For example, you may want to test the i| f 1 1 * 

crashworthiness of cars as a simultaneous function of 4 different design parameters, f | f ± 1 
but with a budget of only three expendable cars. (The issue is not whether this is a 

good plan — it isn't — but rather how to make the best of the situation!) 8 £ ^ of 

The idea is to partition each design parameter (dimension) into M segments, so % J, 1 1 =r 

that the whole space is partitioned into M N cells. (You can choose the segments in |f | g| 

each dimension to be equal or unequal, according to taste.) With 4 parameters and 3 | • § 

cars, for example, you end up with 3x3x3x3 = 81 cells. ? | g |- 3 

Next, choose M cells to contain the sample points by the following algorithm: SJ 5 ^ £ 2 
Randomly choose one of the M N cells for the first point Now eliminate all cells 
that agree with this point on any of its parameters (that is, cross out all cells in the 
same row, column, etc.), leaving (M - 1)^ candidates. Randomly choose one of 

these, eliminate new rows and columns, and continue the process until there is only > |^ % 2 

one cell left, which then contains the final sample point. ? » h 

The result of this construction is that each design parameter will have been f 1 f 8 m 

tested in every one of its subranges. If the response of the system under test is s j| | % % 

dominated by one of the design parameters, that parameter will be found with 9 ^jj ~ o 

this sampling technique. On the other hand, if there is an important interaction | o § 5 g 

among different design parameters, then the Latin hypercube gives no particular § ? | § 5 

advantage. Use with care. - = I S ^ 
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7.8 Adaptive and Recursive Monte Carlo 
Methods 

This section 4 discusses more advanced techniques of Monte Carlo integration. As 
examples of the use of these techniques, we include two rather different, fairly sophisticated, 
multidimensional Monte Carlo codes: vegas [1,2], and miser P]. The techniques that we g. 5 J Q S 
discuss ail fall under the general rubric of reduction of variance (§7.6), but are otherwise < & ^ 
quite distinct 8- & 



^ 6p 



Importance Sampling f If i | 

The use of importance sampling was already implicit in equations (7.6.6) and (7.6.7). i <o 8*3. 

We now return to it in a slightly more formal way. Suppose that an integrand / can be written | 1 i<2"8 

as the product of a function h that is almost constant times another, positive, function g. Then i I f 0< 8 

its integral over a multidimensional volume V is 3 ^ 2. § f 

mil 
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ff<W = f (f/9) 9dV = J hgdV (7.8.1) 

In equation (7.6.7) we interpreted equation (7.8.1) as suggesting a change of variable to 

G y the indefinite integral of g. That made gdV a perfect differential. We then proceeded ^ - « *< ■ 

to use the basic theorem of Monte Carlo integration, equation (7.6.1). A more general 8"g-S 3 pj 

interpretation of equation (7.8.1) is that we can integrate / by instead sampling h — not, z$ •§ ?» =5 

however, with uniform probability density dV, but rather with nonuniform density gdV. In i<j>*Q 3 55 

this second interpretation, the first interpretation follows as the special case, where the means > s.Q <| z 

of generating the nonuniform sampling of gdV is via the transformation method, using the 1 o 3 ^ 

indefinite integral G (see §7.2). §"§ soS 

More directly, one can go back and generalize the basic theorem (7.6.1) to the case § g £5 > 
of nonuniform sampling: Suppose that points Xi are chosen within the volume V with a g | <a 5j 

probability density p satisfying S v,-g Z. Q 

8^3 -—■ ' CO 

pdV=l (7.8.2) f|l| § 

If I jo 3 

The generalized fundamental theorem is that the integral of any function / is estimated, using o | 'n 8 o 
N sample points Xi, , . . , xn, by 3 g S»? g 

o « Z a* 



where angle brackets denote arithmetic means over the N points, exactly as in equation p gf § .|- n 
(7.6.2). As in equation (7.6.1), the "plus-or-minus" term is a one standard deviation error S-k© * * 
estimate. Notice that equation (7.6.1) is in fact the special case of equation (7.8.3), with f §■ i % g 



3. co 



p = constant = 1/V. 

What is the best choice for the sampling density pi Intuitively, we have already ^ 

seen that the idea is to make h = f/p as close to constant as possible. We can be more S o<° w 

rigorous by focusing on the numerator inside the square root in equation (7.8.3), which is > § ^ 

the variance per sample point. Both angle brackets are themselves Monte Carlo estimators § ^ » 

of integrals, so we can write ~ 

5 ■ (£> - <;>* - / - [/ H' - / i w - [/H° < 7 - 8 - 41 

We now find the optimal p subject to the constraint equation (7.8.2) by the functional variation 
6 J lldV-y fdV 2 + X Jpdvj (7.8.5) 



